#include "scene.h"
#include <stdio.h>
#include <math.h>


static short r_grid_gen(struct grid_config *gc, struct grid *g);
static short th_grid_gen(struct grid_config *gc, struct grid *g);
static short z_grid_gen(struct grid_config *gc, struct grid *g);

extern short grid_gen(struct grid_config *gc, struct grid *g)
 {
  short e;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'grid_gen()'...\n");
# endif

  e = r_grid_gen(gc, g);

  if(e >= 0) e = th_grid_gen(gc, g);
  if(e >= 0) e = z_grid_gen(gc, g);
# ifdef TRACE_FILE
  if(e >= 0) fprintf(g->trptr, "...'grid_gen()' ends.\n");
# endif
  return e;
 }


static short r_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *optr;
  int i, cnt;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'r_grid_gen()'...\n");
# endif

# ifdef r_UNIF_GRID
    for(i = 0; i <= g->nr+NrG; i++) g->r[i+NrG] = g->dr*i;
    cnt = g->nr+2*NrG+1;
# else
  g->r[NrG] = 0;
  cnt = NrG+1;
  for(i = 0; i < gc->nr_center; i++, cnt++) g->r[cnt] = (i+1)*(dr_CENTER + 0.5*i*gc->r_cd_center);
  for(i = 0; i < gc->nr_mid12; i++, cnt++) g->r[cnt] = gc->r_mid1 + (i+1)*(gc->dr_mid1 + 0.5*i*gc->r_cd_mid12);
  for(i = 0; i < gc->nr_base; i++, cnt++) g->r[cnt] = gc->r_mid2 + (i+1)*(dr_MID2 + 0.5*i*gc->r_cd_base);
  for(i = 0; i < g->nr_gap; i++, cnt++) g->r[cnt] = r_BASE + (i+1)*(gc->dr_base + 0.5*i*gc->r_cd_gap);

/* Ghost cells beyond the cylindrical boundary. */
  for(i = 0; i < NrG; i++, cnt++) g->r[cnt] = r_CYL + (r_CYL - g->r[NrG+g->nr-i-1]);
# endif

  for(i = 1; i <= NrG; i++) g->r[NrG-i] = -g->r[NrG+i];

  if(gc->r_file != (char *)0)
      {
       optr = (FILE *)fopen(gc->r_file, "w");
       fprintf(optr, "%d\t%E\n", 0, g->r[0]);
       for(i = 1; i < cnt; i++) fprintf(optr, "%d\t%E\t%E\n", i, g->r[i], g->r[i]-g->r[i-1]);
       fclose(optr);
      }

# ifdef TRACE_FILE
  fprintf(g->trptr, "...'r_grid_gen()' ends.\n");
# endif

  return 1;
 }


static short th_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *optr;
  int cnt, i;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'th_grid_gen()'...\n");
# endif

# ifdef th_UNIF_GRID
  for(cnt = 0; cnt <= 2*NthG + g->nth; cnt++) g->th[cnt] = (cnt - NthG)*g->dth;

  if(gc->th_file != (char *)0)
      {
       optr = fopen(gc->th_file, "w");
       for(i = 0; i < cnt; i++) fprintf(optr, "%d\t%E\n", i, g->th[i]);
       fclose(optr);
      }
# else
  fprintf(g->lptr, "ERROR:: in 'th_grid_gen()': No definition for non-uniform grid in theta-direction.\n");
  exit(-1);
# endif

# ifdef TRACE_FILE
  fprintf(g->trptr, "...'th_grid_gen()' ends.\n");
# endif

  return 1;
 }




static short z_grid_gen(struct grid_config *gc, struct grid *g)
 {
  FILE *optr;
  int i, cnt;

# ifdef TRACE_FILE
  fprintf(g->trptr, "Entering into 'z_grid_gen()'...\n");
# endif

# ifdef z_UNIF_GRID
  for(cnt = 0; cnt <= g->nz + 2*NzG; cnt++) g->z[cnt] = (cnt - NzG)*g->dz;
# else

  g->z[NzG] = 0;
  cnt = NzG+1;
  for(i = 0; i < nz_UNIF_BOT; i++, cnt++) g->z[cnt] = g->z[NzG] + (i+1)*gc->dz_unif_bot;
  for(i = 0; i < gc->nz_mid1; i++, cnt++) g->z[cnt] = g->z[NzG+nz_UNIF_BOT] + (i+1)*(gc->dz_unif_bot + 0.5*i*gc->z_cd_bot);
  for(i = 0; i < gc->nz_mid12; i++, cnt++) g->z[cnt] = gc->z_mid1 + (i+1)*(dz_MID1 + 0.5*i*gc->z_cd_mid12);
  for(i = 0; i < gc->nz_top; i++, cnt++) g->z[cnt] = gc->z_mid2 + (i+1)*(dz_MID2 + 0.5*i*gc->z_cd_top);

/* Creating ghost-cells at bottom and top. */
  for(i = 0; i < NzG; i++, cnt++) g->z[cnt] = z_TOP + (z_TOP - g->z[NzG+g->nz-1-i]); /* At the top */
  for(i = 1; i <= NzG; i++) g->z[NzG-i] = g->z[NzG]-g->z[NzG+i]; /* At the bottom */

# endif

  if(gc->z_file != (char *)0)
      {
       optr = fopen(gc->z_file, "w");
       for(i = cnt-1; i >= 1; i--) fprintf(optr, "%d\t%E\t%E\n", i-2, g->z[i], g->z[i]-g->z[i-1]);
       fprintf(optr, "%d\t%E\n", -2, g->z[0]);
       fclose(optr);
      }

# ifdef TRACE_FILE
  fprintf(g->trptr, "...'z_grid_gen()' ends.\n");
# endif

  return 1;
 }
